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Abstract In this paper an autonomous analytical system of ordinary differential 
equations is considered. For an asymptotically stable steady state x° of the system 
a gradual approximation of the domain of attraction (DA) is presented in the case 
when the matrix of the linearized system in a; is diagonalizable. This technique is 
based on the gradual extension of the " embryo" of an analytic function of several 
complex variables. The analytic function is the transformed of a Lyapunov func- 
tion whose natural domain of analyticity is the DA and which satisfies a linear 
non-homogeneous partial differential equation. The equation permits to establish 
an "embryo" of the transformed function and a first approximation of DA. The 
"embryo" is used for the determination of a new "embryo" and a new part of the 
DA. In this way, computing new "embryos" and new domains, the DA is grad- 
ually approximated. Numerical examples are given for polynomial systems. For 
systems considered recently in the literature the results are compared with those 
obtained with other methods. 



1 Introduction 



The domain of attraction (DA) of a steady state of a dynamical system is the 
set of initial states from which the system converges to the steady state itself. In 
order to guarantee stable behavior of a dynamical system in a region of the state 
parameters it is important to know the DA [1]. 

Theoretical research shows that the DA and its boundary are complicated sets 
[2], [3], [4], [5]. In most cases, they do not admit an explicit elementary repre- 
sentation. For this reason, different procedures are used for the approximation of 
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the DA with domains having a simpler shape. This practice became fundamen- 
tal in the last 20 years [6]. The domain which approximates the DA is defined 
by a Lyapunov function, generally quadratic. For a given Lyapunov function the 
computation of the optimal estimate of the DA amounts to solving a non-convex 
distance problem [7], [8], [9], [10], [11], [12], [13]. 

In this paper we present a new technique for a gradual approximation of the 
DA in the case of autonomous analytical systems of ordinary differential equations. 
More precisely for the gradual approximation of the DA of an asymptotically stable 
steady state in which the matrix of the linearized system is diagonalizable. This 
technique is based on the theoretical results obtained in [14], [15], [16] and in the 
following we present a summary of this theory. 



2 Theoretical summary 

We consider the system of differential equations 

x = f{x) (1) 
where / : R™ — > W 1 is an analytical function with the following properties: 

i. /(0) = (i.e. zero is a steady state of (Q)) 

ii. the real parts of the eigenvalues of §j(0) are negative (i.e. x — is an 
asymptotically stable steady state) 

The following theorem provides a tool of determining the DA of the zero steady 
state of (JTJ) . 

Theorem 1. (see [14]) The DA of the null solution of (Qj coincides with the 
natural domain of analyticity of the unique solution V of the problem 

<vv,/) = -H 2 (9) 

V(0) = 1 > 

The function V is positive on DA and lim V(x) = oo for any xo £ dDA. 

X — >Xq 

Thus, the problem of finding the DA is reduced to the determination of the 
natural domain of analyticity of the solution V of (||) . This function will be called 
in the followings the optimal Lyapunov function. 

To determine the optimal Lyapunov function V and its domain of analyticity 
is not easy, but, in the diagonalizable case, we can determine the coefficients of the 
expansion of W = V a S in 0, where S reduces f^(0) to the diagonal form. Then, 
using a Cauchy-Hadamard type theorem, we can obtain the domain of convergence 
D° of the series of W, and DA = S(D°) is a part of the domain of attraction. 
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2.1 The coefficients of the transformed optimal Lyapunov 
function in the diagonalizable case 

For the system ([!]) the following theorem holds: 

Theorem 2. (see [15]) For each isomorphism S : C™ — » C" and g = S^ 1 o f o S, 
the problem 

(VW,g)=-\\Sz\\ 2 

W(0) = o [6 > 

has a unique analytical solution, namely W = V o S where V is the optimal Lya- 
punov function for ^). 

The function W will be called the transformed optimal Lyapunov function. 
In the fallowings, we will suppose that the matrix fj(0) is diagonalizable. 
Let be S : C" — > C" one isomorphism which reduces §~(0) to the diagonal 
form S^ 1 §£(0)5 = diag(Xi, \2...X n ) and g — 5 _1 o/oS. 
For this 5, we consider the expansion of W in 0: 



W(z 1 ,z 2 ,...,z n ) ^2 B h) 



m=2 \ j\=m 

and the expansions in of the scalar components <?i of g: 



(4) 



i (z 1 ,z 2 ,...,z n ) = X t y l + 2J ^2 h )i 

m=2 \j\=m 



jn—jn 1 2 " ™ 



(5) 



Theorem 3. (see [16]) The coefficients Bj 1 j 2 ...j n of the development (Qj are given 
by the following relations: 



■d-E4 '/./ J., i 

i=l 

n 

I] s^s^ £/ |j| = 2 and j p = j q = 1 



(6) 



b'l-i 

E E E[(?i-** + i) 

X^j'iAj p=2 |fc|=p,*i<ji i=l 

,-Sjl-fcl...ii-fci+l...j'n-fcn] if\j\ >3 



According to (||), the coefficients of the terms of degree m > 3 are obtained in 
function of the coefficients of the terms of degree smaller than m. 
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2.2 The domain of convergence of the series of W 

We define the domain of convergence of the series (||) as the set of those z £ C™ 
which have the property that the series (Q) is absolutely convergent in a neighbor- 
hood of z [17]. We denote by D° this domain. Actually, the domain of convergence 
D° coincides with the interior of the set of all the points z° in which the series (^) 
is absolutely convergent (see [17]). 

Theorem 4. ( Cauchy-Hadamard see [17]) A point z belongs to D° if and only if 



lim \B hh ... j A4 2 -zk\<l (7) 

y \j\=rn 

Using D° we can find a part of the DA. 

Theorem 5. (see [15]) If z belongs to D° then the point x = Sz is in the domain 
of attraction DA. 

Actually, by the previous theorem, we obtain the subdomain S(D°) C DA, 
which has the property that dS(D°) n 3D A ^ 0. We also observe that the subdo- 
main S(D°) is symmetrical to the origin (actually, the domain of convergence D° 
is symmetrical to the origin and the axes of coordinates). This first approximation 
of the domain of attraction DA will be denoted by DA = S(D°). 

Condition (^) represents an algoritmizable criterion for determining the region 
of convergence D°. 



2.3 A method of extending the estimate of DA 

In practice, we can compute the coefficients Bj 1 j 2 ___j n up to a finite degree p. This 
degree p has to be big enough for assuming that the domain D® given by 

D° p = {z G C n / IJ2 \B hh ... jn zt4 2 -zt\ < 1} 
V b\=p 

approximates the region of convergence D° of the series of W. On this domain D® 
we approximate W by the "embryo" 

w°(zi,z 2 ,...,z n ) = j2 B hh...3n4 % 4'»-*k (8) 

m=2 \ j\=m 

The first estimate of the region of attraction DA will be DA® — S(D®). 

In order to extend the first estimate DA®, we will expand W in a point z° close 
to dD® in which \Wp(z°)\ is still small. That is because, according to Theorem 
1, the points z close to dD® for which |VFp(z)| is extremely high are close to 
dS- x {DA). 

To find the expansion of W in z° close to dD® , we will compute the expansion 
in z° of the "embryo" W° of W. We obtain: 
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=E E ^,i,...i„(-i - - - (9) 

m=0 \j\=m 

We consider the set 



Dl = {ze C"/ /E - z°)ii(z 2 - z 2 °)^...(z„ - z°H < 1} 

V b'l=f 

which provides a new part DAp = S(D p ) of DA 

So, DA^ 1 U DA*, gives a larger estimate of the DA We can continue this 
procedure for a few steps, till the values \W p \ become extremely large and we 
obtain the estimate DA° p U DA\ U ... U DA k p = S(D° p U D* U ... U D£) of DA 

3 Numerical results 

The computations were made using a program written in Mathematica 4, Wolfram 
Research. The data for the estimations (the degree up to which the approximation 
is made, the necessary timing for the estimations) are displayed in Table 1. 

3.1 Systems with known domains of attraction 

In this subsection, we will present some examples of systems of two or three dif- 
ferential equations, for which we can compute easily the DA We will apply our 
technique to these examples, and we will show how the real domains of attraction 
are gradually approximated. These examples are meant to validate our procedure. 

3.1.1 Example 1 

This is an example of a system for which the null solution has a bounded domain 
of attraction: 




The domain of attraction of the null solution of this system is the interior of the 
circle of radius 2 centered in (1, 0): 

DA = {(x 1 ,x 2 ) e R 2 /( Xl - if +x 2 2 < 4} 

After three steps, we obtain the estimate shown in Figure 1. 
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The thick black line represents the true boundary of the domain of attraction, 
the dark grey set denotes the first estimate DA® of DA and the further estimates 
DAp of DA with k > 1 are colored in light grey. 




Figure 1: The estimate of DA obtained after three steps for system ( J10| ) 
3.1.2 Example 2 

The following system of three differential equations is considered: 

X\ = — Xi(l — x\ — x\ + £§) 

x\ = -X 2 (l - x\ - x\ +xf) (11) 

x 3 = -a; 3 (l - x\ - x\ + xl) 
The boundary of the DA of the null solution of this system is: 

3D A = {(xi,x 2 ,x 3 ) e M 3 /l - x\ - x% + x\ = 0} 

The first estimate DA® is shown in Figure 2.1. After 2 steps, we obtain the 
estimate shown in Figure 2.2. 




Figure 2.1: The estimate of DA obtained Figure 2.2: The estimate of DA obtained 
after 1 step for system ( |TT| ) after 2 steps for system ( pT| ) 
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3.2 Systems for which we don't know the DA 

In this subsection, some systems of differential equations are presented for which 
we don't know the DA. For these examples, we will apply the technique presented 
in Section 2, and we will show that the estimate obtained using this technique is 
better than the estimates obtained in [10]. 



3.2.1 Example 3 

In [10], the following example is considered: 

{ fl = X2 (12) 
1 x\ = —2xi — 3x2 + x\ — x\ + x\Xi 

In Figure 3.1, an estimate of the DA is shown obtained after two steps for three 
different points close to the boundary of DA® . We observe that this estimate 
covers the estimate presented in [10]. Figure 3.2 presents the estimate of DA 
obtained after four steps. The thick black lines plotted in the following figures 
represent a part of the approximated boundary of the DA. The black interrupted 
line represents the boundary of the estimate of the domain of attraction obtained 
in [10]. 




Figure 3.1: The estimate of DA after 2 
steps for three different points close to 
the boundary of DA® for system (|12j) 



Figure 3.2: The estimate of DA obtained 
after 4 steps for system ( [l2"|) 



3.2.2 Example 4 

In [10], the following system of three equations is considered: 



(13) 

- 3x 2 - 2x 3 + x\ + x\ + X3 

Figure 4.1 shows the estimate of DA obtained after one step. The estimate of the 
DA obtained after two steps is presented in Fig. 4.2. 
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Figure 4.1: The estimate of DA obtained Figure 4.2: The estimate of DA obtained 
after 1 step for system ( |l3| ) after 2 steps for system ( p"3| ) 



Table 1. Numerical data 



example 


order of approximation 


timing for I s * step 


timing for 2 nd step 


1 


50 


1.7 min 


34.7 min 


2 


30 


10.2 min 


14.7 h 


3 


30 


0.9 min 


9.9 min 


4 


30 


19.1 min 


16.2 h 
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